! SK SOLID WALL BCS 
subroutine usdef3( un, u,  ub, ucb, ucsb  &
                 , pn,p, pb,pcb, pcsb &
                 , ron,rho,rhob,rhocb,rhocsb &
                 ,ds1,esm,dt,dtmin,epsilon2, NoSlipCoef)
	!IMPLICIT REAL *8(a-h,o-z)
use Data_cfoam, only: epsMaximumPrinciple
implicit none

    real(8),intent(in) :: NoSlipCoef
    real *8 esm(3)
	real *8  un(3), u(3),  ub(3), uf(3), ucb(3), ucsb(3), ucf(3), ucsf(3)
	real *8 pn,p, pb, pf, pcb, pcsb, pcf, pcsf
	real *8 rhon,rho,rhob,rhof,rhocb,rhocsb,rhocf,rhocsf
	real *8 ds1,ds2,dt,sound,rofon,dtmin
	real *8 ut(3,9)
	integer j
    real *8 ei,z1,vtmp,dttmp
    real *8 char3,q,r,qn,rn,u1x,ron,vf,wf,vb,wb
    real *8 char3f,soundf,gf,sf,qf,qcf,qcsf,char2f,qwell,qmaxf,qminf
    real *8 char3cf
    real *8 char3csf
    real *8 char3b,soundb,gb,sb,rb,rcb,rcsb,char2b,rwell,rmaxb,rminb
    real *8 char3cb
    real *8 char3csb
    real *8 charc1,charc2,charc3
    real *8 a(3,3)
    real *8 epsilon
    real *8 es(3)
    real *8 defb
    real *8 char3in,dtmint,gam,gam2,gamx,gamy,gin,pin,qin,roin,s,scb,scsb,sn, &
        sin,uin,vin,win,cv
    external defb
    real *8 epsilon2
    
    real(8) :: alfa, sn1,sn2, ut12,ut13, ut22,ut23

    epsilon=epsilon2
    !if (pcsb.lt.0.95d5) epsilon=0.2
    !if (pcsb.gt.1.15d5) epsilon=0.2
    gam=1.4d0
    gamx=2.0d0*dsqrt(gam)/(gam-1.0d0) ! A
    gamy=((gam-1.0d0)/(2.0d0*gam))
    gam2=2.0d0*gam
    es=esm/dsqrt(esm(1)**2+esm(2)**2+esm(3)**2)
	call coor(es(1),es(2),es(3),a)
! 1-u 2-ub 3-uf 4-ucb 5-ucsb 6-ucf 7-ucsf
    call  transform(u,ut(1,1),a)
    call  transform(ub,ut(1,2),a)
 
    call  transform(ucb,ut(1,4),a)
    call  transform(ucsb,ut(1,5),a)
 

    char3 =ut(1,1)
    char3b=ut(1,2)
 
    char3cb=ut(1,4) 
    char3csb=ut(1,5)

    s=dlog(p)-gam*dlog(rho)
    sb=dlog(pb)-gam*dlog(rhob)

    scb=dlog(pcb)-gam*dlog(rhocb)
    scsb=dlog(pcsb)-gam*dlog(rhocsb)


	gb=gamx*dexp(scsb/gam2)
	gf=gb

    soundb=dsqrt(gam*pcsb/rhocsb)

	r=char3+gb*p**gamy
	rb=char3b+gb*pb**gamy
	rcb=char3cb+gb*pcb**gamy
	rcsb=char3csb+gb*pcsb**gamy
	char2b=char3csb+soundb
	rwell=(rcsb-rcb)/(0.5d0*dt)+char2b*(r-rb)/ds1
	rmaxb=dmax1(r,rb,rcb)+dt*rwell + epsMaximumPrinciple
	rminb=dmin1(r,rb,rcb)+dt*rwell - epsMaximumPrinciple

    rn=(2d0*rcsb-rb*(1.0d0-epsilon))/(1.0d0+epsilon)
	if(rn>rmaxb) rn=rmaxb
	if(rn<rminb) rn=rminb
! SK DEFINE THE INVARIANTS THROUGH THE CHARACTERISTIC BCS
	qn=-rn
	u1x=0
!	(rn+qn)/2.
	pn=((rn)/gb)**(2.0d0*gam/(gam-1.0d0))
! 	pn=((rn-qn)/(gb+gf))**(2.0d0*gam/(gam-1.0d0))

    call GetAlfa(char3csb,alfa)
	if(alfa>=1.0d0)then ! inflow to boundary face
	    sn=scsb

	    ut(1,8)=0.0d0
	    ut(2,8)=ut(2,5)
	    ut(3,8)=ut(3,5)
	else if(alfa<=0.0d0)then ! ouflow from boundary face
	    sn=defb(s,sb,scb,scsb,char3csb,dt,ds1,epsilon)
	    ut(1,8)=0.0
	    ut(2,8)=defb(ut(2,1),ut(2,2),ut(2,4),ut(2,5),char3csb,dt,ds1,epsilon)
	    ut(3,8)=defb(ut(3,1),ut(3,2),ut(3,4),ut(3,5),char3csb,dt,ds1,epsilon)
	else
	    sn1=scsb
	    ut12=ut(2,5)
	    ut13=ut(3,5)

	    sn2=defb(s,sb,scb,scsb,char3csb,dt,ds1,epsilon)
	    ut22=defb(ut(2,1),ut(2,2),ut(2,4),ut(2,5),char3csb,dt,ds1,epsilon)
	    ut23=defb(ut(3,1),ut(3,2),ut(3,4),ut(3,5),char3csb,dt,ds1,epsilon)
	    
	    sn = alfa * sn1 + (1.0d0 - alfa) * sn2
	    ut(1,8) = 0.0d0
	    ut(2,8) = alfa * ut12 + (1.0d0 - alfa) * ut22
	    ut(3,8) = alfa * ut13 + (1.0d0 - alfa) * ut23
	end if

    ut(1,8) = NoSlipCoef * ut(1,8)
    ut(2,8) = NoSlipCoef * ut(2,8)
    ut(3,8) = NoSlipCoef * ut(3,8)

!    un(1) = 0.0
!    un(2) = 0.0
!    un(3) = 0.0

    call transformr(ut(1:3,8),ut(1:3,9),a)
    do j=1,3    
        un(j)=ut(j,9) ! +u1x*es(j) ! u1x = 0.0
    end do

!	if(char3csb.le.0) then ! inflow to boundary face
!	    sn=scsb
!
!	    un(1)=ucsb(1)-char3csb*es(1)
!	    un(2)=ucsb(2)-char3csb*es(2)
!	    un(3)=ucsb(3)-char3csb*es(3)
!
!!	    ut(1,8)=0.0
!!	    ut(2,8)=ut(2,5)
!!	    ut(3,8)=ut(3,5)
!!       call transformr(ut(1:3,8),ut(1:3,9),a)
!
!	else ! ouflow from boundary face
!	    sn=defb(s,sb,scb,scsb,char3csb,dt,ds1,epsilon)
!	    ut(1,8)=0.0
!	    ut(2,8)=defb(ut(2,1),ut(2,2),ut(2,4),ut(2,5),char3csb,dt,ds1,epsilon)
!	    ut(3,8)=defb(ut(3,1),ut(3,2),ut(3,4),ut(3,5),char3csb,dt,ds1,epsilon)
!        
!        call transformr(ut(1:3,8),ut(1:3,9),a)
!	    do j=1,3    
!            un(j)=ut(j,9)+u1x*es(j)
!	    end do
!	end if
	
	ron=(pn/exp(sn))**(1d0/gam)
	vtmp=dsqrt(un(1)**2+un(2)**2+un(3)**2)
    sound=sqrt(gam*pn/ron)
    dtmint=ds1/(sound+vtmp)
	if (dtmint.lt.dtmin) dtmin=dtmint

    return 
end subroutine usdef3